import pandas as pd
import numpy as np
from itertools import product
import matplotlib.pyplot as plt
import seaborn as sns
%matplotlib inline
sns.set_context('notebook')
from IPython.display import display
drop_pos = ['dropoff_longitude','dropoff_latitude']
pick_pos = ['pickup_longitude', 'pickup_latitude']
lon_pos = ['pickup_longitude', 'dropoff_longitude']
lat_pos = ['pickup_latitude', 'dropoff_latitude']
position = ['pickup_longitude', 'pickup_latitude','dropoff_longitude','dropoff_latitude']
week = ['Sunday', 'Monday', 'Tuesday', 'Wednesday','Thursday', 'Friday', 'Saturday']
train = pd.read_csv('train.csv')
test = pd.read_csv('test.csv')
display(train.head())
display(train.describe().style)
A partir da descrição do pandas, podemos notar, observando os campos min, max, e os quartis 1 e 3:
def add_date_specifics(df, column_name, new_prefix):
df.loc[:,column_name] = pd.to_datetime(df[column_name])
df[new_prefix +'year'] = df[column_name].dt.year
df[new_prefix +'month'] = df[column_name].dt.month
df[new_prefix +'yearday'] = df[column_name].dt.dayofyear
df[new_prefix +'hour'] = df[column_name].dt.hour
df[new_prefix +'minute'] = df[column_name].dt.minute
df[new_prefix +'weekday'] = df[column_name].dt.weekday_name
df[new_prefix +'weeknum'] = df[column_name].dt.weekday
df[new_prefix +'weekend'] = train['pickup_weekday'].isin(['Saturday', 'Sunday'])
def add_distances(df):
"""
Abaixo usamos a aproximação euclideana, ignorando a curvatura da terra,
por tratarmos de uma unica cidade
"""
ny_lat = 40.7 # From wikipedia
R = 6371 # From wikipedia
deg_rad_ratio = np.pi/180
conversion = np.array([R*deg_rad_ratio, R*deg_rad_ratio*np.cos(ny_lat*deg_rad_ratio)])
df['line_distance'] = np.sqrt((((df[pick_pos].values - df[drop_pos].values)*conversion)**2).sum(axis=1))
df['manh_distance'] = np.abs(((df[pick_pos].values - df[drop_pos].values)*conversion)).sum(axis=1)
#
#df.apply(lambda x: great_circle(x[pick_pos], x[drop_pos]), axis=1)
def add_speed(df):
""" Velocidade em km/h"""
df['speed(km/h)'] = 3600*df['line_distance']/df['trip_duration']
df['manhattan speed(km/h)'] = 3600*df['manh_distance']/df['trip_duration']
whole = train.append(test)
add_date_specifics(train, 'pickup_datetime', 'pickup_')
add_date_specifics(train, 'dropoff_datetime', 'dropoff_')
add_distances(train)
add_speed(train)
train['log_trip_duration'] = np.log10(train['trip_duration']+1)
add_date_specifics(test, 'pickup_datetime', 'pickup_')
add_distances(test)
train['duration_from_date'] = (train['dropoff_datetime']-train['pickup_datetime']).dt.total_seconds().astype(int)
print('Number of duration inconsistencies:',
(train.duration_from_date != train.trip_duration).sum(),
'(time provided vs time from datetime diff)')
display('Train:',train.isnull().sum(), 'Test: ', test.isnull().sum())
Podemos ver que não há valores faltando nem no conjunto de treino nem no de teste.
plt.figure(figsize=(12,4))
plt.subplot(1,3,1)
sns.distplot(np.log10(train['line_distance']+1), kde=False); plt.xlabel('log (line_distance)');
plt.title('Distancia das viagens em log')
plt.subplot(1,3,2)
plt.title('Distancia das viagens curtas < 1km')
sns.distplot(train['line_distance'][train['line_distance']<1], kde=False);
plt.subplot(1,3,3)
plt.title('Distancia das viagens curtas < .1km')
sns.distplot(train['line_distance'][train['line_distance']<.1], kde=False);
plt.tight_layout()
Notamos que existem muitas viagens que não mudam de lugar, provavelmente problemas na utilização do dispositivo pelos taxistas.
Removeremos as viagens com tamanho < .1km, que é aproximadamente 1 quarteirão. As viagens de tamanho 0 serão removidas juntamente.
train = train[train['line_distance']>.1]
plt.figure(figsize=(12,4))
plt.subplot(1,3,1)
sns.distplot(np.log10(train['speed(km/h)']+1)); plt.xlabel('log( speed(km/h) )');
plt.subplot(1,3,2)
ax=sns.distplot(train['speed(km/h)'][train['speed(km/h)'] > 200], kde=False);
plt.subplot(1,3,3)
ax=sns.distplot(train['speed(km/h)'][train['speed(km/h)'] < 200], kde=False);
De acordo com esse site e esse, a velocidade maxima em Nova York é 65mph (~104km/h), somente nas rodovias rurais.
Achamos razoável assumir um limite máximo de 200km/h para o nosso dataset, levando em conta que as distâncias medidas são aproximadas. Essa decisão é bastante conservadora, visto que a distância em linha reta tenderá a subestimar as distâncias e as velocidades, e 200km/h é uma velocidade bastante difícil de ser alcançada.
train = train[train['speed(km/h)'] < 200]
train['log_line_distance'] = np.log10(train['line_distance'])
train['log( speed(km/h) )'] = np.log10(train['speed(km/h)'])
plt.figure(figsize=(18,6))
ax=plt.subplot(1,4,1)
sns.barplot(x='passenger_count', y=0, hue='vendor_id', data=train.groupby(['vendor_id','passenger_count']).apply(len).reset_index(), ax = ax, palette = 'Pastel2')
ax=plt.subplot(1,4,2)
sns.violinplot(x='passenger_count', y='log_line_distance', hue='vendor_id', data=train, ax = ax, palette = 'Pastel2')
ax=plt.subplot(1,4,3)
sns.violinplot(x='passenger_count', y='log_trip_duration', hue='vendor_id', data=train, ax = ax, palette = 'Pastel2')
ax=plt.subplot(1,4,4)
sns.violinplot(x='passenger_count', y='log( speed(km/h) )', hue='vendor_id', data=train, ax = ax, palette = 'Pastel2')
plt.tight_layout()
Podemos ver que há poucas viagens com 0 passageiros, e que as distribuições são semelhantes aos outros tipos de viagem.
As viagens com passageiros > 6 foram removidas em passos anteriores, da distância mínima e velocidade.
plt.figure(figsize=(12,10))
plt.subplot(3,3,1)
sns.distplot(train['log_trip_duration']); plt.title('Duração das viagens')
plt.subplot(3,3,2)
sns.distplot(train['trip_duration'][train['trip_duration']<60], kde=False); plt.title('Duração das viagens < 1min (zoom)')
plt.subplot(3,3,3)
sns.distplot(train['log_trip_duration'][train['log_trip_duration']>4.4], kde=False);
plt.title('Duração das viagens > 7h (zoom)'); plt.ylim(0,50);
plt.subplot(3,3,4)
sns.distplot(np.log10(train['speed(km/h)']+1)); plt.xlabel('log( speed(km/h) )');
plt.subplot(3,3,5)
sns.distplot( train['speed(km/h)'][train['trip_duration'] < 20], kde=False );
plt.title('Velocidade nas viagens curtas')
plt.subplot(3,3,6)
ax=sns.distplot( train['speed(km/h)'][train['log_trip_duration']>4.4], kde=False )
plt.title('Velocidade nas viagens longas')
#plt.xlabel('log( speed(km/h) )');
plt.subplot(3,3,7)
sns.distplot( train['log_line_distance'], kde=True )
plt.subplot(3,3,8)
sns.distplot( train['line_distance'][train['trip_duration'] < 20], kde=False )
plt.title('Distancia das viagens curtas')
plt.subplot(3,3,9)
sns.distplot( train['line_distance'][train['log_trip_duration']> 4.4], kde=False )
plt.title('Distancia das viagens longas')
plt.tight_layout()
Não parece razoável haverem tantas viagens com duração próxima a 1 dia, (log_trip_duration = 5), mas a duração > 10 dias é completamente não razoável, e provavelmente houve problema no dispositivo, que ficou ligado sem parar.
train = train[train.log_trip_duration<6]
plt.figure(figsize=(12,4))
plt.subplot(1,3,1)
sns.boxplot(x='variable', y='value', data=train[lon_pos].melt())
plt.subplot(1,3,2)
sns.boxplot(x='variable', y='value', data=train[lon_pos].melt()); plt.ylim(-85, -60)
ax=plt.subplot(2,3,3)
train[lon_pos][train[lon_pos]>-90].quantile(np.arange(100)/100).plot(ax=ax)
ax=plt.subplot(2,3,6)
train[lon_pos][train[lon_pos]>-90].quantile(np.arange(20)/100).plot(ax=ax,legend=False);
ax=ax.twinx()
train[lon_pos][train[lon_pos]>-90].quantile(np.arange(10,100)/100).plot(ax=ax,legend=False);
plt.tight_layout()
plt.figure(figsize=(32,16))
plt.subplot(1,2,1)
plt.scatter(train['pickup_longitude'], train['pickup_latitude'], s=.3);
plt.xlim(-74.05, -73.77); plt.ylim(40.55,40.93); plt.title('Pickup')
plt.subplot(1,2,2)
plt.scatter(train['dropoff_longitude'], train['dropoff_latitude'], s=.3);
plt.xlim(-74.05, -73.77); plt.ylim(40.55,40.93); plt.title('Dropoff')
Manualmente e por iterações, foram decididos limites para as coordenadas de latitude e longitude que mantivessem as caracteristicas visualmente importantes no mapa, e que serão usadas nos mapas a partir daqui.
No mapa dos pickups é fácil notar dois 'bumps' fora da cidade de Manhattan. Uma pesquisa rápida no google maps mostra que se tratam do aeroporto regional de La Guardia, mais próximo à cidade, e do aeroporto internacional Jonh F. Kennedy, mais distante.
min_lon, max_lon = -74.05, -73.77
min_lat, max_lat = 40.55, 40.93
min_pos = pd.Series(dict(zip(position,[min_lon,min_lat, min_lon, min_lat])))
max_pos = pd.Series(dict(zip(position,[max_lon,max_lat, max_lon, max_lat])))
is_central = ((train[position] > min_pos) & (train[position] < max_pos)).sum(axis=1)==4
plt.figure(figsize=(10,10))
ax = plt.subplot(2,2,1)
sns.distplot(train[is_central]['pickup_longitude'])
ax = plt.subplot(2,2,2)
sns.distplot(train[is_central]['dropoff_longitude'])
ax = plt.subplot(2,2,3)
sns.distplot(train[is_central]['pickup_latitude'])
ax = plt.subplot(2,2,4)
sns.distplot(train[is_central]['dropoff_latitude'])
plt.tight_layout()
quants = pd.DataFrame()
for p in np.arange(0,1.001,.001):
quants = quants.append(train[position].quantile(p))
quants[lon_pos].reset_index().set_index(['pickup_longitude'])['index'].plot(label='pickup longitude')
quants[lon_pos].reset_index().set_index(['dropoff_longitude'])['index'].plot(linestyle='--',label='dropoff longitude')
plt.legend(loc='center right')
plt.twiny()
quants[lat_pos].reset_index().set_index(['pickup_latitude'])['index'].plot(color='g',label='pickup latitude')
quants[lat_pos].reset_index().set_index(['dropoff_latitude'])['index'].plot(linestyle='--',color='k', label='dropoff latitude')
plt.ylabel('Cumulative probability')
plt.legend(loc='upper left')
plt.title('')
É possível ver que ainda possuimos outliers muito distantes tanto para latitude quanto longitude.
plt.figure(figsize=(18,8))
ax=plt.subplot(2,3,1)
sns.barplot(x='vendor_id', y=0, data= train.groupby('vendor_id').apply(len).reset_index(), palette = 'Pastel2')
plt.ylabel('occurence')
ax=plt.subplot(2,3,2)
sns.violinplot(x='passenger_count', y='log_trip_duration', hue='vendor_id', data=train, ax = ax, palette = 'Pastel2')
ax=plt.subplot(2,3,3)
sns.barplot(x='passenger_count', y=0, hue='vendor_id', data=train.groupby(['vendor_id','passenger_count']).apply(len).reset_index(), ax = ax, palette = 'Pastel2')
ax=plt.subplot2grid((2,3),(1,0))
sns.barplot(x='pickup_hour', y=0, hue='vendor_id', palette = 'Pastel2',
data = train.groupby(['vendor_id', 'pickup_hour']).apply(len).reset_index())
plt.title('Number of drives per hour of day')
ax=plt.subplot2grid((2,3),(1,1))
sns.barplot(x='pickup_hour', y='trip_duration', data=train, ax=ax, hue='vendor_id', palette = 'Pastel2')
plt.title('Mean duration of trips')
ax=plt.subplot2grid((2,3),(1,2))
sns.barplot(x='pickup_hour', y='speed(km/h)', data=train, ax=ax, hue='vendor_id', palette = 'Pastel2')
plt.title('Mean speed of trips')
plt.tight_layout()
plt.figure(figsize=(18,10))
# Weekday-------------
ax=plt.subplot2grid((3,3),(0,0))
train.pickup_weekday.value_counts()[week].plot.bar(ax=ax, rot=15)
plt.title('Number of drives per day of week'); plt.ylim(150000); plt.ylabel('Number of drives')
ax=plt.subplot2grid((3,3),(1,0))
sns.barplot(x='pickup_weekday', y='trip_duration', data=train, order=week, ax=ax)
plt.title('Mean duration of trips'); plt.ylim(400);
ax=plt.subplot2grid((3,3),(2,0))
sns.barplot(x='pickup_weekday', y='speed(km/h)', data=train, order=week, ax=ax)
plt.title('Mean speed of trips');
# Hour -------------
ax=plt.subplot2grid((3,3),(0,1))
sns.barplot(x='pickup_hour', y=0, hue='pickup_weekend', palette = 'Set3',
data = train.groupby(['pickup_weekend', 'pickup_hour']).apply(len).reset_index())
plt.title('Number of drives per hour of day'); plt.ylabel('Number of drives')
ax=plt.subplot2grid((3,3),(1,1))
sns.barplot(x='pickup_hour', y='trip_duration', data=train, ax=ax, hue='pickup_weekend', palette = 'Set3')
plt.title('Mean duration of trips'); plt.ylim(400);
ax=plt.subplot2grid((3,3),(2,1))
sns.barplot(x='pickup_hour', y='speed(km/h)', data=train, ax=ax, hue='pickup_weekend', palette = 'Set3')
plt.title('Mean speed of trips');
# Month-------------
ax=plt.subplot2grid((3,3),(0,2))
sns.barplot(x='pickup_month', y=0,
data=train.groupby('pickup_month').apply(len).reset_index(), palette='Set2')
plt.title('Number of drives per month of year'); plt.ylim(150000); plt.ylabel('Number of drives')
ax = plt.subplot2grid((3,3),(1,2))
sns.barplot(x='pickup_month', y='trip_duration', data=train, ax=ax, palette='Set2')
plt.title('Mean duration of trips'); plt.ylim(400);
ax=plt.subplot2grid((3,3),(2,2))
sns.barplot(x='pickup_month', y='speed(km/h)', data=train, ax=ax, palette='Set2')
plt.title('Mean speed of trips');
plt.tight_layout()
É possivel notar uma diferença clara entre os dias da semana e os fins de semana, entre as 5 da manha e as 18. As viagens costumam ser muito mais lentas nos dias da semana, o que provavelmente se deve ao commute, da ida e volta diária ao trabalho.
Sem atenção às horas, ainda é possivel ver uma diferença de velocidade nas viagens, mais lentas no meio da semana, e que curiosamente se mantém altas nas segundas feiras, mesmo sendo um dia de trabalho.
def citymap(df, pos = pick_pos, func=len, precision=4000):
lats = np.arange( int(precision*train[is_central]['dropoff_latitude'].min()),
int(precision*train[is_central]['dropoff_latitude'].max())+1)
lons = np.arange( int(precision*train[is_central]['dropoff_longitude'].min()),
int(precision*train[is_central]['dropoff_longitude'].max())+1)
densemap = np.zeros((lons.shape[0], lats.shape[0]))
local = df.copy()
local[pos] = local[pos].apply(np.around, decimals=4)
sparse_vals = local.groupby(pos).apply(func)
for (lon, lat), count in sparse_vals.iteritems():
i = int(precision*lon-lons[0])
j = int(precision*lat-lats[0])
if i < densemap.shape[0] and j < densemap.shape[1] and i>=0 and j>=0:
densemap[i,j] = count
return pd.DataFrame(densemap, index = pd.Index(lons/precision, name='Longitude'), columns=pd.Index(lats/precision, name='Latitude'))
pickmap = np.log10(citymap(train[is_central]) +1)
dropmap = np.log10(citymap(train[is_central], pos = drop_pos) +1)
plt.figure(figsize=(20,30))
ax = plt.subplot(2,1,1)
sns.heatmap(pickmap.T[::-1]); plt.title('Pick-ups')
ax = plt.subplot(2,1,2)
sns.heatmap(dropmap.T[::-1]); plt.title('Drop-offs');
del(pickmap)
del(dropmap)
Observando os mapas de pickup e dropoff, e possivel notar que existe muito mais consistencia nos pickups, o que e bastante razoavel, levando em conta que os taxistas devem pegar passageiros em vias movimentadas e em seus respectivos pontos de taxi, e provavelmente nao passam muito tempo procurando passageiros em vias pequenas. Essa analise considera que o meio mais comum de pegar taxi, para essas duas empresas em questao, ainda e "presencial", em oposiçao a atraves de aplicativos.
Outra possibilidade e que os passageiros tenham menor tendencia a pegar taxis na ida para a cidade, talvez pela maior facilidade do transporte publico nessa direçao, e peguem mais comumente na volta para casa, ou em viagens relacionadas a trabalho dentro da propria cidade. Essa hipotese bate com o fato de que o numero de taxis e muito maior nas horas entre 7h e 23h, em uma diferença que ocorre principalmente nos dias da semana, enquanto no fim de semana a quantidade de taxis e menos variavel.
O mais provavel e que ambos os fatores exerçam alguma influencia sobre os dados.
plt.figure(figsize=(80,12))
for i in range(8):
plt.subplot(1,8,i+1)
to_plot = (train['dropoff_hour']==i+5) & is_central
plt.scatter(train[to_plot]['dropoff_longitude'], train[to_plot]['dropoff_latitude'],
s=5, alpha=.4, c=np.log10(train[to_plot]['speed(km/h)']), cmap='plasma', vmin=.7, vmax=1.7)
plt.title('Dropoff at %dam'%(i+5))
plt.axis('off')
plt.figure(figsize=(80,12))
for i in range(8):
plt.subplot(1,8,i+1)
to_plot = (train['dropoff_hour']==i+16) & is_central
plt.scatter(train[to_plot]['dropoff_longitude'], train[to_plot]['dropoff_latitude'],
s=5, alpha=.4, c=np.log10(train[to_plot]['speed(km/h)']), cmap='plasma', vmin=.7, vmax=1.7)
plt.title('Dropoff at %dpm'%(i+4))
plt.axis('off')
É possível perceber que, enquanto a densidade de viagens terminando no centro de Manhattan diminue com a chegada da noite, ela aumenta perceptivelmente nos arredores, indicando provavelmente a volta de trabalhadores para suas casas.
A velocidade média das viagens, indicada pela cor dos pontos, indica uma melhora gradual no transito, mais rapida na periferia, especialmente mais distante da cidade entre os aeroportos. Às 23h ainda sao lentas as viagens para a parte mais central da cidade.
raw_feat_list = ['vendor_id','pickup_hour','pickup_month','pickup_weeknum', 'store_and_fwd_flag',
'pickup_weekend', 'line_distance', 'manh_distance'] + position
onehot_feat_list = ['pickup_hour', 'pickup_month', 'pickup_weekday']
from sklearn.preprocessing import OneHotEncoder
onehot = OneHotEncoder(sparse=False)
def one_hot_df(df, feature):
onehot = OneHotEncoder(sparse=False)
dense_X = onehot.fit_transform(df[feature].values.reshape(-1,1))
return pd.DataFrame( dense_X, columns = ['%s_%d'%(feature, i) for i in onehot.active_features_] )
X = pd.get_dummies(train[onehot_feat_list], columns = onehot_feat_list).join(train[raw_feat_list])
X['store_and_fwd_flag'] = X['store_and_fwd_flag'] == 'Y'
y = train.log_trip_duration
Xtest = pd.get_dummies(test[onehot_feat_list], columns = onehot_feat_list).join(test[raw_feat_list])
Xtest['store_and_fwd_flag'] = Xtest['store_and_fwd_flag'] == 'Y'
from xgboost import XGBRegressor, plot_importance
from sklearn.model_selection import cross_val_score
def rmsle(y_true,y_pred):
"""
From https://www.kaggle.com/wiki/RootMeanSquaredLogarithmicError
"""
assert len(y_true) == len(y_pred)
return np.square(np.log(y_pred + 1) - np.log(y_true + 1)).mean() ** 0.5
def rmse(y_true,y_pred):
assert len(y_true) == len(y_pred)
return np.square(y_pred - y_true).mean() ** 0.5
from sklearn.metrics import make_scorer
rmsle_score = make_scorer(rmsle)
rmse_score = make_scorer(rmse)
from skopt import gp_minimize
def optimize(X, y, n_calls=80, n_random_starts=10):
hyperparameter_space = [(1,15),
(0.01,1.),
(10,300),
(1e-6,10.,"log-uniform"),
(1,30),
(1e-6,1.,"log-uniform"),
(1e-6,1.,"log-uniform"),
(1e-2,1.,"log-uniform"),
(1e-2,1.,"log-uniform")]
hyperparameter_names = ['max_depth', 'learning_rate', 'n_estimators',
'gamma', 'min_child_weight', 'reg_alpha', 'reg_lambda',
'colsample_bytree', 'subsample']
def objective_(params):
parameters = dict(zip(hyperparameter_names, params))
print(parameters)
score = np.mean(np.log(10)*cross_val_score(XGBRegressor(**parameters), X, y, cv=3, scoring=rmse_score))
return score
return gp_minimize(objective_, hyperparameter_space, n_calls=n_calls,n_random_starts=n_random_starts, n_jobs=3,verbose=3)
%%time
res = optimize(X, y)
best = {'max_depth': 11, 'learning_rate': 0.21476942062372031, 'n_estimators': 287,
'gamma': 2.3074723022921413e-06, 'min_child_weight': 16, 'reg_alpha': 0.5815787682811846,
'reg_lambda': 0.60169547700624582, 'colsample_bytree': 0.78028984074973284, 'subsample': 0.36043975312396048}
xgb = XGBRegressor(**best)
xgb.fit(X, y)
ypred = xgb.predict(Xtest)
fig, ax = plt.subplots(1, figsize=(15,15))
plot_importance(xgb, ax=ax)

sub = pd.DataFrame(test.id)
sub['trip_duration'] = 10**ypred
sub.to_csv('my_sub.csv',index=False)

Ao fim das análises e submissão do kaggle, notei shortcomings que deveriam ser corrigidos para melhoria do modelo:
Os citado acima seriam os previstos próximos passos na análise desse dataset.
As maiores dificuldades encontradas nesse trabalho foram: